solo: --soloUMIfiltering MultiGeneUMI_All (closes #144) - #152

CONTRIBUTING.md requires the description to match the code; the entries for the
other themes split out of scverse#152 belong to their own PRs.

CONTRIBUTING.md requires the description to match the code; the entries for the
other themes split out of scverse#152 belong to their own PRs.

* fix: --runThreadN 1 ran on every core, not on one

The rayon global pool was configured only when `--runThreadN` was
greater than 1. Skipping the build at 1 does not give one thread: it
leaves rayon's default, which is one worker per logical core. So
`--runThreadN 1` ran the whole machine.

Measured on 200k reads, before: 2.44 s wall at **1300% CPU**. After:
26.23 s wall at 100% CPU. The old figure was not a fast single-threaded
run, it was a sixteen-way run wearing the wrong flag.

This matters beyond the flag reading falsely. A scheduler or a container
given one CPU gets sixteen worker threads; on a shared machine the run
oversubscribes every other job; and with a thread-caching allocator each
of those threads keeps its own heap, which is the very cost the comment
above this code says the pool sizing exists to avoid.

It also means the project's thread-invariance checks were weaker than
they read: the `--runThreadN 1` leg was not a one-thread leg. Verified
now that it is one: records are byte-identical between 1 and 8 threads
on 200k real reads, and byte-identical to the previous binary's output
at `--runThreadN 1`. Only the `@PG` `CL:` line differs between thread
counts, because it records the command line.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(changelog): record the --runThreadN 1 fix

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* refactor: drop Transcript::read_seq, which nothing ever reads

The field was filled with a full copy of the read at every finalised
transcript and assigned again in four more places, and no code in the
crate reads it. Not "reads it rarely": the compiler was asked, and after
deleting the field every one of the 97 resulting errors is a struct
literal or an assignment. There is no read site, in `src` or in tests.

Measured on 200k real reads at 8 threads: about 50 transcripts are
finalised per read, so removing it takes **15.1 million allocations off a
200 million total** (7.6%) and 3.87 GB of copying.

Wall clock does not move, and that is worth recording rather than hiding:
six interleaved rounds at 87-92% CPU idle give medians 20.90 s against
20.65 s with the direction mixed, inside the run-to-run spread. mimalloc
is fast enough that seventy-five small allocations per read do not
surface. The reason to remove it is that it is dead weight, not that it
is slow.

That number also calibrates #168 downward: if removing 7.6% of the
allocations changes nothing measurable, the rest of the allocation
programme is unlikely to be worth a new dependency.

`Transcript::read_seq` is `pub`, so this is an API removal and needs
sign-off. Nothing outside the crate can be relying on its contents being
meaningful, though, since it is only ever written.

Output-neutral: SAM byte-identical on 200k real reads.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(changelog): record the Transcript::read_seq removal

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* test: an A/B harness that refuses to average over a busy machine

Every perf measurement on this project needs the same three guards, and
I got each of them wrong at least once in a single session:

An earlier harness sampled CPU idle once, before the first run. Three
separate measurements then drifted below the threshold mid-series and it
kept going, leaving me to spot the contaminated rounds by eye in the
output. This samples idle before and after every run, drops a round if
any of the four samples falls short, and prints how many it dropped, so
a median over four surviving rounds cannot be mistaken for a median over
six.

The check is on CPU idle rather than load average, because load average
is an exponential average over minutes: it refused to measure at 2.24 on
a machine whose cores were all free.

It reports the spread within each side next to the difference between
the medians, and says so in as many words when the difference is
smaller. Two changes I measured looked like wins on medians alone and
were inside the spread.

Both sides run as ./rustar-aligner with --outFileNamePrefix ./ from
inside their own directory, because the @pg CL: line records argv
verbatim: running ./old against ./new is enough to make the output
differ, which cost me two false "output is not neutral" alarms.

The header documents the fourth trap, which no script can enforce:
timing a total hides the part that changed. BAM writing is 1-4% of a
yeast run, so a total dominated by alignment cannot resolve a change to
the writer. Run the None configuration alongside and read the
difference.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* feat(solo): --soloCellReadStats CB writes CellReads.stats

One row per cell barcode with the fifteen counters STARsolo reports: how the
barcode matched, whether the read mapped to one locus or several, whether it
landed on a feature, where in the gene and on which strand, whether it was
mitochondrial, and whether it reached the matrix. The per-cell UMI and gene
totals come from the raw matrix rather than from the read counters, so they
agree with what the matrix says by construction.

Reads whose barcode never resolved are summed into a single `CBnotInPasslist`
row instead of being dropped. That row is the reason the file is useful: it is
the difference between "these cells look thin" and "most of the input never
reached a cell at all".

The region columns split by strand — an antisense read counts under `exonicAS`
or `intronicAS`, never under `exonic` or `intronic`.

`--genomeChrSetMitochondrial` names the chromosomes behind the `mito` column.
Without it the column is zero throughout, which is honest: no chromosome was
declared mitochondrial.

D24 comes with it. STAR emits these rows by walking a libc++ `unordered_map`,
which at these sizes is the reverse of each barcode's first appearance. That is
reproduced, including across threads: the per-read accumulator merges in read
order, so a threaded run writes the same file as a serial one. It stops being
reproducible past the point where libc++ rehashes, since the order then depends
on the bucket count. The values never differ, only which line they sit on.
Recorded in docs-old/dev/divergences.md.

* docs: record the CellReads.stats row order in DIVERGENCE.md

Section 3.2, in the format CONTRIBUTING.md asks for.

* docs(changelog): keep only this PR's entry

CONTRIBUTING.md requires the description to match the code; the entries for the
other themes split out of #152 belong to their own PRs.

* refactor(solo): drop CellReadStats::merge, which nothing calls

Reads are folded in under a mutex, so there are no per-thread partials to merge;
the function was reachable only from its own test. CONTRIBUTING.md rules out
shipping a function no production path reaches, and the PR description claimed
its test as evidence of thread-safety that the mutex actually provides.

* feat(solo): --runMode soloCellFiltering

Cell-calls an existing raw count matrix without aligning anything, taking the
raw directory and an output prefix as STAR does:
`--runMode soloCellFiltering /path/to/raw/ /path/to/out/prefix`.

Cell calling is a decision about a matrix, not about reads. Re-calling with
different `--soloCellFilter` parameters should not mean re-aligning 400 million
reads, and a matrix produced by another tool should be callable too.

The matrix is streamed into the same temp-body form the align path builds, so
`called_cells` and `emptydrops_called` are the identical code here and there
rather than a second implementation free to drift. Counts are rounded on the way
in: a multimapper matrix carries real values, and the filters work on UMI
totals.

`--runMode` becomes a token list, because that is what STAR's is: the mode
followed by its arguments. The mode itself is now validated rather than falling
back to `alignReads`, so a typo is refused instead of quietly running something
else.

The standalone `emptydrops` binary still exists and still carries its own copy
of the algorithm, which no longer matches this one. Removing it means moving
`test/solo_genefull_compare.py` and `test/solo_genefull_h5_compare.py` to the
new mode first, so it is left alone here rather than broken.

---------

Co-authored-by: Benjamin Demaille <benjamin.demaille@icloud.com>
Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>

* fix(solo): MultiGeneUMI_CR gives a tied UMI to nobody, not to everybody

`--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest
read count. CellRanger's rule is the opposite on exactly that case: the
gene with the *strictly* highest count takes the UMI, and a tie means no
gene counts it.

STAR walks the genes keeping a running maximum and clears its winner
whenever it meets an equal count
(`SoloFeature_collapseUMIall.cpp:212-224`):

    if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; }
    else if (ig.second==maxu) { maxg=-1; };
    ...
    if ( maxg+1==0 ) continue; // not counted for any gene

One read per gene is the ordinary shape of a multi-gene UMI, and it is
always a tie, so the old rule made the flag inert in practice rather
than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells
from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared
between two genes), against STAR 2.7.11b with the same flags:

                       identical entries    STAR counts   rustar counts
    before             13 749 / 14 806          15 423          16 465
    after              13 902 / 13 967          15 423          15 414

The flag removed nothing at all before; STAR removes 1 030 counts. The
gap goes from +1 042 to -9.

The outcome does not depend on the order the genes are visited — a
strict maximum always ends as the winner, a tie always ends with none —
so iterating a `HashMap` here stays deterministic.

`multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests
missed: they only covered 3 reads against 1, where both rules agree.

Not yet implemented, and stated so rather than left to be discovered:
STAR applies a second condition, that the winning gene must also hold
the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file,
lines 226-232). That needs the pre-correction counts, which this code
does not keep. The 65 entries still differing out of 13 967 are the
place to look for its effect.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(changelog): record the MultiGeneUMI_CR tie fix

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* fix(solo): MultiGeneUMI_CR decides ownership on corrected UMIs

STAR corrects UMIs within each gene *before* deciding which gene owns a
UMI, and applies two conditions, not one
(`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`):

1. one gene must hold a strictly higher read count than every other, on
   the **corrected** UMI map — that is #173, already landed;
2. and that winner must not be beaten in the **uncorrected** map at the
   same key.

The second condition exists because correction moves reads between UMIs:
a gene can win only because correction folded a neighbouring UMI onto it,
and STAR rejects that win rather than counting it.

Reproducing it needs the order STAR uses. The generic path here filters
multi-gene UMIs first and corrects afterwards, which cannot express either
condition: by the time correction happens the ownership decision is
already made. `MultiGeneUMI_CR` therefore takes its own path, which is
also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`,
so there is no combination this bypasses.

`cellranger_1mm_map` exposes the correction mapping that
`cellranger_1mm` already computed and threw away.

Measured against **CellRanger 10.0.0** on the 20 000-read fixture from
#172, with #165 and #173 also applied:

                              identical entries   CellRanger   rustar
    #165 + #173               13 651 / 13 709        15 111     15 091
    plus this change          13 676 / 13 709        15 111     15 116

Entries CellRanger has and we do not go from 29 to 7, and the count gap
from -20 to +5, which is 0.03%.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(divergence): drop a reference to a test that no longer exists

* feat(solo): bit-exact libc++ mt19937, generate_canonical and discrete_distribution

STARsolo's `EmptyDrops_CR` rescue draws from `std::mt19937`, converts to
doubles with `std::generate_canonical<double, 53>`, and picks categories with
`std::discrete_distribution`. Two of those three are implementation-defined in
the parts that matter: the standard fixes mt19937's output but not how
`generate_canonical` consumes it, and says nothing about how
`discrete_distribution` maps a uniform onto categories.

So porting "the algorithm" is not enough — it has to be libc++'s algorithm,
because that is what STAR is built against and where its numbers come from.
libc++ accumulates two 32-bit draws in *ascending* significance and divides by
2^64; a most-significant-first accumulation, or one draw scaled to 53 bits,
both give perfectly good uniforms and neither reproduces STAR.

Every expected value in the tests came out of a C++ program compiled against
the real libc++ and run, not from reading its source. `tests/libcxx_oracle.cpp`
is that program, kept so the values can be regenerated rather than trusted.
`generate_canonical` is compared as bit patterns, since a difference in the
last place changes which category a sample lands in.

Not yet wired into the EmptyDrops path. `solo::count` samples with a
`SplitMix64` stream under a comment calling it "WeightedIndex-equivalent;
empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in
general, since two unrelated generators cannot agree on an arbitrary number of
draws. It is true of whatever was checked and unknown elsewhere. Replacing it
moves cell calls, so it belongs in its own change with the solo differential
run against it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* feat(solo): EmptyDrops_CR uses Simple Good-Turing and libc++'s sampler

Two approximations in the CellRanger cell-calling path are replaced by what
CellRanger and STAR actually compute. Both move cell calls, which is the point:
the previous numbers were plausible rather than right.

The ambient profile is now smoothed with Simple Good-Turing (Gadsby & Sampson,
via Elworthy's implementation, which is what STAR vendors). The ambient counts
come from a small sample of empty droplets, so a gene seen twice there is not
twice as likely as one seen once, and a gene seen zero times is not impossible —
it is one the sample was too small to show. SGT fits the frequency spectrum and
reserves mass for the unseen from the singleton rate, then smooths the rest
along a log-log line. What was here before had the right shape and the wrong
numbers: it reserved mass the same way but distributed the remainder in
proportion to raw counts, with no smoothing at all.

The Monte-Carlo null is now drawn with libc++'s `std::mt19937` and
`std::discrete_distribution`, seeded `19760110 * (isim + 1)` per simulation, as
STAR seeds it. The previous sampler was a SplitMix64 stream under a comment
calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell
calls" — a claim that cannot hold in general, since two unrelated generators
cannot agree over an arbitrary number of draws. The libc++ types were ported and
checked against real libc++ in the previous commit on this branch; this wires
them in. One generator per simulation, no shared state, so the walks still run
in any order on any number of threads and give the same p-values.

D17 comes with it: STAR leaves `PZero` uninitialised when the spectrum has fewer
than five distinct frequencies and `analyse()` bails, so it reads whatever the
stack held. Here it is zero from construction, which is what "no basis for
reserving unseen mass" means. Recorded in docs-old/dev/divergences.md.

* docs: record the EmptyDrops SGT divergence in DIVERGENCE.md

Section 1.2, in the What STAR does / What rustar-aligner does / Why / Impact /
Source format CONTRIBUTING.md asks for, replacing the docs-old file the earlier
version of this work carried.

* docs(divergence): file the EmptyDrops entry under section 1, note the
second RNG

* fix(params): refuse MultiGeneUMI_CR without --soloUMIdedup 1MM_CR

* feat(solo): --soloFeatures Transcript3p, with --soloClusterCBfile

Quantifies transcripts rather than genes, from where each read's 3' end sits
relative to each transcript's. In a 3'-biased assay that distance is what
separates isoforms: a read 200 bases from the end of one and 4000 from the end
of another is evidence for the first. The distribution of those distances is
estimated from the run's own histogram, smoothed and cut where the 3' peak
decays into the body, and used as the likelihood in an EM over UMIs.

Concordance needed no new code. `align_to_transcripts` already refuses to
project an alignment that leaves the transcript, touches an intron, or crosses a
junction the transcript does not have — which is exactly STAR's `Concordant`
(`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant;
one that does not, is not. The projection also puts the 5' end at coordinate
zero for both strands, so the distance to the 3' end is one expression rather
than two.

Two behaviours worth stating because they are not the obvious ones:

Output is per cluster, not per cell, and `--soloClusterCBfile` is required.
A single cell does not have enough UMIs to resolve isoforms, so the EM would be
fitting noise. Asking for the feature without a clustering is refused rather
than run.

A UMI seen on several reads contributes the *intersection* of their transcript
sets. Those reads came from one molecule, so a transcript missing from any of
them cannot be its source. Taking the union would let a single stray read
resurrect an isoform every other read excluded.

Two of STAR's quirks are reproduced rather than corrected, because the cut point
and every weight depend on them: the running-average divisor is `min(2N+1,
i + N)` rather than the number of elements actually summed, and the transcript
length factor is taken from the cumulative distribution at `trLen - 1`
(`SoloFeature_quantTranscript.cpp`).

Numbers are formatted the way C++'s default stream prints them — six
significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since
the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s
default disagree on both notation and digit count.

* refactor(solo): drop Transcript3pAcc::merge, which nothing calls

Records are accumulated under a mutex, so there are no partials to merge. It was
dead from the moment it was written; CONTRIBUTING.md rules out shipping it.

* docs(solo): note that STAR marks Transcript3p under development

parametersDefault puts both Transcript3p and --soloClusterCBfile between
"#####UnderDevelopment_begin : not supported - do not use" and
"#####UnderDevelopment_end", and STAR --help prints that banner around
them. The module said none of this.

It matters for how the port is read: it follows STAR's code, so it
inherits the unfinished parts of that code, and a differential against
STAR compares two implementations of something STAR does not support.
A reviewer should be told that before deciding to take it.

* fix(solo): implement MultiGeneUMI_All instead of aliasing it to MultiGeneUMI

`--soloUMIfiltering MultiGeneUMI_All` resolved to the same variant as
`MultiGeneUMI`, which is neither what STAR does nor what the option is
documented to do. Of the three available behaviours it was the only one nobody
had asked for.

In STAR the option is a no-op: it is parsed and stored, but its consumption site
tests only the `MultiGeneUMI` flag, so selecting it leaves the filter entirely
off. Documented, it removes a UMI seen in more than one gene from *all* of them,
rather than from the losers only.

`UmiFiltering::MultiGeneUmiAll` now exists and does the documented thing: a UMI
appearing in several genes is evidence of a collision or of chimeric
amplification, so it is discarded outright rather than attributed to whichever
gene happened to read deepest. Single-gene UMIs are untouched, which the test
checks across every mode.

Raised upstream as #144 before changing it, since "be faithful to STAR" and "do
what the flag says" genuinely point in opposite directions here.

Also adds `docs-old/dev/divergences.md`, recording this and the homopolymer-UMI
rule, so deliberate differences are written down rather than rediscovered as
surprises in a differential run.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(changelog): record the MultiGeneUMI_All fix

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

* docs(divergence): correct the MultiGeneUMI_All entry, defer the
homopolymer one

---------

Co-authored-by: Benjamin Demaille <benjamin.demaille@icloud.com>
Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>